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Large-eddy simulations of the turbulent flow in a lid-driven cubical cavity have been carried out at a Reynolds 
number of 12 000 using spectral element methods. Two distinct subgrid-scales models, namely a dynamic 
Smagorinsky model and a dynamic mixed model, have been both implemented and used to perform long-lasting 
simulations required by the relevant time scales of the flow. All filtering levels make use of explicit filters applied 
in the physical space (on an element-by-element approach) and spectral (modal) spaces. The two subgrid-scales 
models are validated and compared to available experimental and numerical reference results, showing very 
good agreement. Specific features of lid-driven cavity flow in the turbulent regime, such as inhomogeneity 
of turbulence, turbulence production near the downstream corner eddy, small-scales localization and helical 
properties are investigated and discussed in the large-eddy simulation framework. Time histories of quantities 
such as the total energy, total turbulent kinetic energy or helicity exhibit different evolutions but only after a 
relatively long transient period. However, the average values remain extremely close. 
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I. INTRODUCTION 

The study of a lid-driven flow of a Newtonian fluid in a 
rectangular three-dimensional cavity is of particular interest 
in view not only of the simplicity of the flow geometry but 
also the richness of the fluid flow physics manifested by multi- 
ple counter-rotating recirculating regions at the corners of the 
cavity depending on the Reynolds number, Taylor-Gortler-like 
(TGL) vortices, flow bifurcations and transition to turbulence. 
This flow structure is now well documented thanks to a rela- 
tively rich literature reporting both computational and experi- 
mental studies. A comprehensive review of the fluid mechan- 
ics of driven cavities is provided by Shankar and Deshpande 
infill. 

In the present paper, our focus resides in relatively high- 
Reynolds-number and three-dimensional lid-driven cubical 
cavity flows. At Reynolds number higher than a critical value 
comprised between 2 000 and 3 000, an instability appears in 
the vicinity of the downstream corner eddy B BE- As the 
Reynolds number further increases, turbulence develops near 
the cavity walls, and at Reynolds number higher than 10 000, 
the flow near the downstream corner eddy becomes fully tur- 
bulent. The highest Reynolds number attained was 12000 
by direct numerical simulation (DNS) performed by Leriche 
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and Gavrilakis JU and 10000 experimentally by Koseff and 
Street |0] and Prasad and Koseff 0]. In the literature, papers 
using the lid-driven cavity problem as a benchmark test case 
to evaluate the performance of numerical algorithms are pro- 
liferating, but are often limited to two space dimensions or 
to Reynolds numbers below 10 000. More recently, one may 
however notice the important developments of novel and more 
physical numerical methods applied to the lid-driven cavity 
flow such as molecular dynamics by Chen and Lin in ^ and 
also the lattice-Boltzmann model applied by He et al. in jfj. 

The results reported herein correspond to the numerical 
simulation of the flow in a lid-driven cubical cavity at the 
Reynolds number of 12 000 placing us in the locally-turbulent 
regime. The spatial discretization relies on spectral element 
methods (SEM) which have been mainly applied to the DNS 
of fluid flow problems at low and moderate Reynolds num- 
bers. With the advent of more powerful computers, espe- 
cially through cluster technology, higher Re values seem to 
fall within the realm of feasibility. However, despite their high 
accuracy, spectral element methods are still far from reaching 
industrial applications that involve developed turbulence at Re 
values of the order of 10 6 — 10 7 . The reason for that dismal 
performance is that a resolved DNS including all scales from 
the largest structures to Kolmogorov scales, needs a number 
of degrees of freedom that grows like Re 9 / 4 . Therefore with 
increasing Re, one has to increase the number of elements, E, 
and the degree, N, of the polynomial spaces. This places the 
computational load far out of the reach of present day comput- 
ers. Large-eddy simulation (LES) represents an alternative to 
DNS insofar that it involves less degrees of freedom because 
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the behavior of the small scales are modeled. The numerical 
simulations presented in this paper encompass two different 
LES based on two distinct subgrid-scales modeling both us- 
ing an eddy-viscosity assumption, and one using in addition a 
mixed model relying on the scale-similarity hypothesis, simi- 
larly to Zang et al. in fiHEl for Re = 10 000. Compared to 
the previous works of Zang et al, the two LES reported here 
offer simulation length ten times larger therefore increasing 
the accuracy of the ensemble averaging and more importantly 
allowing to capture intermittent turbulent production. These 
events lead to the determination of large eddies suggested to 
be mainly responsible for the turbulence production near the 
downstream corner eddy. 

Unlike low-order methods such as finite volumes or finite 
differences, spectral and spectral element methods allow a 
complete decoupling between the mathematical formulation, 
the subgrid modeling, the numerical technique and the filter- 
ing technique, which are introduced successively in Sec. [TT] 
Specifically, we are first seeking to validate the two subgrid- 
scale models introduced in Sec. |II]which rely on explicit filter- 
ing techniques specific to spectral element spatial discretiza- 
tion. Sec. Hn] presents a short, but comprehensive valida- 
tion procedure. In Sec. II VI emphasis is put on characterizing 
the turbulent flow in its locally-turbulent regime. Fundamen- 
tal features are qualitatively and quantitatively investigated 
such as the inhomogeneity of the turbulence, the turbulence 
production in the downstream-corner-eddy region, the small- 
scales turbulent structures in the cavity flow and finally the 
peculiar helical properties. 



II. THE MODEL AND NUMERICAL TECHNIQUE 
A. Mathematical modeling 




FIG. 1 : Sketch of the geometry of the lid-driven cubical cavity. 

The fluid enclosed in the cavity is assumed to be incom- 
pressible, Newtonian with uniform density and temperature. 
The flow is governed by the Navier-Stokes equations inside 
the fluid domain denoted by V = [—h, +h] 3 with no-slip 
boundary conditions on every cavity walls, except on the top, 



see Fig. Q] The flow is driven by imposing a prescribed veloc- 
ity distribution with nonzero mean on the "top" wall — named 
lid in the sequel — with the velocity field maintained every- 
where parallel to a given direction. The details regarding the 
imposition of this Dirichlet boundary condition for the veloc- 
ity field at the lid is discussed in Sec. Ill C 31 As the flow 
presents turbulent zones coexisting with laminar regions, the 
numerical simulation incorporates the mathematical models 
involved by the large-eddy simulation method in order to re- 
solve the complex dynamics of the flow. As a consequence, 
the governing equations of the large-eddy simulations are the 
filtered Navier-Stokes equations. Large-scale quantities, des- 
ignated in the sequel by an "overbar", are obtained by a filter- 
ing procedure on the computational domain V = [—1, +1] 3 
using h for the non-dimensionalization of lengths. The appli- 
cation of a low-pass inhomogeneous and anisotropic spatial 
filter to the Navier-Stokes equations in the Eulerian velocity- 
pressure formulation and in convective form for the nonlinear 
term yields 

d\\ 

— +U- Vu = -Vp + z/Au- V ■ r, (1) 

at 

V ■ u = 0, (2) 

in which u is the filtered velocity field, t denotes the time, 
p = P/p is the filtered static pressure and v the assuredly 
constant, uniform kinematic viscosity. The symbols V and A 
represent the nabla and Laplacian operators, respectively. The 
subgrid-scale (SGS) stress tensor r is given by 

t = uu — uu, (3) 

and accounts for the effects of the unresolved- or small-scales 
on the dynamics of the resolved- or large-scales lfl2tl . 

B. Subgrid-scale models 

1. Under-resolved direct numerical simulation 

In the same framework as it prevails among the praction- 
ers, one can resort to the DNS computations without any 
LES model, but with the nodal filtering technique described 
in Sec. HID ll to let the numerical method dissipate locally the 
high-wave-number modes introduced by the insufficient space 
discretization. Such an approach corresponds to an under- 
resolved DNS (UDNS). 

2. Smagorinsky model 

The SGS Smagorinsky model (SM) lfl3tl uses the concept 
of turbulent viscosity and assumes that the small scales are 
in equilibrium, balancing energy production and dissipation. 
This yields the following expression for the eddy viscosity 

^ sgs = (C s A) 2 |S|, (4) 

where |S| = (2E ij S ij ) 1 / 2 is the magnitude of the filtered 
strain-rate-tensor with Sij = l/2(diii/dxj + duj/dxt), Cs 
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is the Smagorinsky constant and A the filter width. The 
Smagorinsky model has several drawbacks. The most severe 
one is the constant value of Cs during the computation which 
produces too much dissipation. Furthermore the SM does not 
provide the modeler with backscattering where kinetic energy 
is transferred from small scales to larger scales in an inverse- 
cascading process. 



3. Dynamic Smagorinsky model 

The dynamic Smagorinsky model (DSM) proposed by Ger- 
mano et al. 03] overcomes the difficulty of constant Cs, by 
allowing it to become dependent of space and time. Now we 
have a dynamic parameter Cd = Cd(x, t). Let us introduce 
a test-filter length scale A that is larger than the grid length 
scale A (e.g. A = 2A). Using the information provided by 
those two filters and assuming that in the inertial range of the 
turbulence energy spectrum, the statistical self-similarity ap- 
plies, we can better determine the features of the SGS stress. 
As detailed in with the test filter, the former LES Eq. © 
yields a dynamic parameter having the expression 



C d 



where 



(a -13): L d 
(a - 3) : (a - 3) ' 



L = uu uu, 



a = -2A |S|S, 







-2A ISIS. 



(5) 

(6) 

(7) 
(8) 



The notation ":" in Eq. (0 is used for inner tensor product 
(double contraction), and the upper index d denotes the devi- 
atoric part of the tensor. 



4. Dynamic mixed model 



The dynamic mixed model 111 111 introduced to tackle cav- 
ity flows is a blend of the mixed model of Bardina et al. 
IU6I1 and the former dynamic Smagorinsky model. We notice 
that Bardina's scale similarity model is not an eddy-viscosity 
based model. Instead it belongs to the class of structural mod- 
els lfl2ll and relies on the scale-similarity principle. It pro- 
duces almost no dissipation and for that reason needs to be 
used jointly with dissipative models such as the Smagorin- 
sky model — Bardina's mixed model — or with the dynamic 
Smagorinsky model. The approach of Zang et al. IU1I1 was 
extended by Liu et al. lfl7ll who proposed a new similarity 
subgrid-scale model for incompressible flows, in which the 
subgrid stress tensor is assumed to be proportional to the re- 
solved stress tensor. Vreman et al. 111811 later modified the 
DMM formulation to remove a mathematical inconsistency by 
expressing the scale-similarity part of the sub-test-scale stress 
T — see lfl5ll — using only u. Salvetti and Banerjee lfl9ll and 



Horiuti (20tl extended the DMM to two distinct dynamic two- 
parameter models. Morinishi and Vasilyev ll2"lll recommended 
a modification to the dynamic two-parameter mixed model of 
Salvetti and Banerjee lfl9fl for large-eddy simulation of wall 
bounded turbulent flow. The works of Vreman et al. ll22ll and 
Winckelmans et al. l23ll also closely relate to the DMM ap- 
proach. As mentioned by Morinishi and Vasilyev ll2lll and 
Ghosal [24], the reliability of the results of large-eddy simu- 
lation is strongly affected by both the effectiveness of the sub- 
grid scale model and the accuracy of the numerical method, 
particularly in the approximation of the non-linear convective 
term. As mentioned in Sec. Q] the SEM is decoupled from 
the subgrid modeling and offers a high accuracy characteris- 
tics of spectral methods. Therefore, the present work focuses 
on the one-parameter type of dynamic mixed model DMM as 
introduced by Zang et al. Il ill for the lid-driven cavity flow. 
The modification suggested by Vreman et al. lfl8ll was not 
implemented; a priori tests with their modified DMM using 
samples from the DNS results by Leriche and Gavrilakis [5] 
showed no noticeable improvement over the DMM of Zang 
et al. in the subgrid stress correlations. Therefore increasing 
the computational expense by adding an additional filtering 
level operation as required by the modification of Vreman et 
q/. lflUl . seemed unjustified. 

By decomposing the velocity field as 

u = u + u', (9) 

where u' represents the subgrid-scale velocity field and by in- 
serting in Eq. J3j, we can redefine the SGS stress as proposed 
by Germano 12511 



t = C + C + K, 



where 



C = 



u u — u u, 



C = u u' + u' u — (u u' + u' u) , 



(10) 



(ID 



K = u'u' u' u', 



are designated as the modified Leonard stress, the SGS cross 
term, and the modified SGS Reynolds stress, respectively. The 
modified Leonard term can be calculated by resolved quanti- 
ties and corresponds essentially to the mixed model. The two 
other terms are unresolved residual stresses and are treated 
through the Smagorinsky model, see ifTHl for greater details. 
Following the same dynamic procedure as in lfl5ll and Sec. 
Ill B 31 and with the same notations, one obtains a dynamic 
coefficient which reads 



C d = 



where 



(a - 3) : (G d - L d ) 
(a - 3) : (a - 3) ; 



u u — u u. 



(12) 



(13) 



The expression of the dynamic coefficient given in Eq. dT2b 
for the dynamic mixed model is similar to the one for the dy- 
namic model — see Eq. (f5]) — the tensor L d being replaced by 

g d - L d . 
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C. Numerical technique 



discrete Navier-Stokes equations become 



1. Space discretization 

The numerical method treats Eqs. (HJ-© within the weak 
Galerkin formulation framework. The spatial discretization 
uses Lagrange-Legendre polynomial interpolants. The reader 
is referred to the monograph by Deville et al. l26ll for full 
details. The velocity and pressure are expressed in the Pjy — 
Pjv_2 functional spaces where Pjv is the set of polynomials of 
degree lower than N in each space direction. This spectral el- 
ement method avoids the presence of spurious pressure modes 
as it was proved by Maday and Patera lf27l l28ll . The quadra- 
ture rules are based on a Gauss-Lobatto-Legendre (GLL) grid 
for the velocity nodes and a Gauss-Legendre grid (GL) for the 
pressure nodes. 

Borrowing the notation from l26ll . the semi-discrete filtered 
Navier-Stokes equations resulting from space discretization 
are 



,du 



M— + Cu + i/Ku — D p + Dt = 0, (14) 
at 



= 0. 



(15) 



The diagonal mass matrix M is composed of three blocks, 
namely the mass matrices M. The global vector u contains 
all the nodal velocity components while p is made of all nodal 
pressures. The matrices K, D T , D are the discrete Laplacian, 
gradient and divergence operators, respectively. The matrix 
operator C represents the action of the non-linear term written 
in convective form u • V, on the velocity field and depends 
on u itself. The semi-discrete equations constitute a set of 
non-linear ordinary differential equations ( fT4b subject to the 
incompressibility condition (TT3T >- 



2. Time integration 

The state-of-the-art time integrators in spectral methods 
handle the viscous linear term and the pressure implicitly by a 
backward differentiation formula of order 2 to avoid stability 
restrictions such that 



uAt < C/N 4 



(16) 



while all non-linearities are computed explicitly, e.g. by a 
second order extrapolation method, under the CFL restriction 
u max At < C/N' 2 . Nonetheless, as the LES viscosity is not 
constant, we modify the standard time scheme in such a way 
that this space varying viscosity be handled explicitly as this 
was done e.g. in lfl5l I29I l3(ill . Let us define the effective 
viscosity as 



ftff 



+ 0, 



eff 



0, 



(17) 



where i/ cst is the sum of the physical viscosity v and the aver- 
age of ^ sgs over the computational domain. The filtered semi- 



M^= + %Ku - D T £ = -Cu + 2D(^ eff - u cst )S, (18) 



Du = 0, 



(19) 



and the previous time splitting still applies. The viscous ex- 
plicit term on the right-hand side does not harm stability as 
the magnitude of the term 2D(i/ e ff — I'csOS is less than that of 
Cu. 

The implicit part is solved by a generalized block LU de- 
composition with a pressure correction algorithm l26l l3lll . 



3. The lid-filtered velocity distribution 

As already mentioned by Leriche and Gavrilakis in 0], im- 
posing a given velocity distribution on the lid of a cavity is 
neither an easy task experimentally nor numerically. Indeed 
imposing a constant lid velocity profile leads to a singularity 
(discontinuous behavior in the velocity boundary conditions) 
at the edges and at the corners of the lid, see Fig. Q] With- 
out adequate treatment, this discontinuous behavior will un- 
dermine the convergence and the accuracy of any numerical 
method in the vicinity of the lid. For the two-dimensional 
case, a well known solution (but with no physical relevance) 
is to subtract the most singular terms of the analytical expres- 
sion of the local stream-function expansion near the lid cor- 
ners. The extension of such procedure to three-dimensional 
cases is still missing even though several recent attempts are 
reported, see I2I I32I1 . In order to explicitly filter the discontin- 
uous behavior, the constant lid velocity profile is regularized 
by the use of a high-order polynomial expansion which van- 
ishes along its first derivatives at the lid edges and corners 



1 2 r 

u(x,y = h, z,t) = Uq l—(x/h) l—(z/h) 



v(x, y = h, z,t) = w(x, y = h, z, t) = 0. 



(20) 



This profile flattens very quickly near the lid edges and cor- 
ners while away from them, it grows rapidly to a constant 
value over a short distance. The exact form and the polyno- 
mial order of the profile is discussed in (0, B|. The highest 
polynomial order of this distribution in both x- and z-direction 
is 36. Such high-order polynomial expansions lead to steep 
velocity gradients in the vicinity of the edges of the lid. The 
grid refinement, in terms of spectral element distribution near 
the lid will be presented in greater details in Sec. HID One of 
the constraint in the grid design is to ensure the proper reso- 
lution of the lid velocity distribution by the spectral element 
decomposition. 



D. Filtering techniques 

As spectral elements offer high spatial accuracy, we con- 
struct explicitly the filters using two spectral techniques. The 
first one is a nodal filter acting in physical space on the nodal 
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velocity components (and pressure) to render the computa- 
tions stable in the long range integration. The second method 
is designed as a modal filter and is carried out in spectral space 
in an element by element fashion. That filter corresponds 
specifically to the convolution kernel of the low-pass LES fil- 
tering. 

1. Nodal filter 

The nodal filter due to Fischer and Mullen ll33ll is ade- 
quately suited to parallel spectral element computation. In- 
troducing h]yj,j = 0, . . . , N the set of Lagrange-Legendre 
interpolant polynomials of degree N based on the GLL grid 
nodes £jv,fe, k = 0, . . . , N, the rectangular matrix operator 
Iff of size (M + 1) x (N + 1) is such that 

{^N)ij = hN,j(^M,i)- (21) 

Therefore, the matrix operator of order N — 1 

ITjv-i = In-i ^at" 1 ' (22) 

interpolates on the GLL grid of degree N — 1 a function de- 
fined on the GLL grid of degree N and transfers the data back 
to the original grid. This process eliminates the highest modes 
of the polynomial representation. The one-dimensional (ID) 
filter is given by the relation 

u = [crfljv-i + (1 - a)I%]u. (23) 

The LES version of the filter sets a = 1 and is given by 

u = iZ ijfu, (24) 

where M is equal to N — 2 or TV — 3. The three-dimensional 
(3D) extension results easily from the matrix tensor product 
properties of the filter. It is worth noting that by construction 
such nodal filter constitutes a projective filter, i.e u = u. 

2. Modal filter 

Here, the variable u is approximated by a modal basis first 
proposed in the p-version of the finite elements and used by 
Boyd rt34ll as a filter technique. It is built up on the reference 
parent element as 

<po = — y~ j 91 = ( 25 ) 

fa = L k (0 - L k _ 2 (®, 2<k<N. 

Conversely to the Lagrange-Legendre nodal basis used in our 
spectral element calculations, this modal basis d25l l forms a 
hierarchical set of polynomials allowing to define in an ex- 
plicit and straightforward manner a low-pass filtering proce- 
dure. The one-to-one correspondence between the nodal La- 
grangian basis and the p-representation yields 

N 

= ^2^k<Pk((,i), (26) 



which in matrix notation reads 

u = *u. (27) 

The low -pass filtering operation is performed in spectral space 
through a diagonal matrix T with components such that 

T =T 1 = 1 and T k = {1 + { l /kc)2) 2 < k < N, 

(28) 

where the cut-off value k c corresponds to T k = 1/2. The 
entire filtering process for a one-dimensional problem is given 
by 

u= G*u = x u. (29) 

The three-dimensional extension is again trivial by the matrix 
tensor product properties. 

As noted in 12911 . the effect of such modal filter onto a given 
field expanded in the modal basis d25l l presents the interesting 
feature of maintaining the inter-element C° -continuity. More 
rigorously, such C°-continuity is enforced if and only if both 
<f>c and <f>i are not at all affected by the low-pass filtering, in 
other words if and only if To = T\ = 1. Nevertheless, it 
has been observed that such C°-breakage does not constitute 
a major issue for our simulations as it only affects the eddy 
viscosity field and other terms present into the modeling of the 
effect of the SGS tensor (01, and which are not used directly 
for constructing a solution retaining the C°-continuity feature. 

Such modal filter is invertible and consequently is not pro- 
jective, i.e u/u. 

3. The filter length 

The decomposition of the computational domain into spec- 
tral elements of given sizes, within which a GLL distribution 
of grid points based on the polynomial degree is chosen, re- 
quires a specific definition of the filter length A. In order to 
account for both the size of each spectral element and its value 
of the polynomial order, and following l35ll . the filter length 
for a ID spectral element method is chosen as 

A = -, (30) 
V 

where s is the element size and p the highest polynomial de- 
gree in the spectral decomposition Eq. (|26| | that is the closest 
to the cut-off frequency k c . In the particular context of the 
modal filter previously introduced, p is such that 

p = k, such that inf(T fe ) < T kc = 1/2, k = 0, . . . , N. 

(31) 

We notice that the filter length decreases when the element is 
refined. The straightforward extension of Eq. d30b to our 3D 
problem using rectilinear elements leads to 

A(x,y,z) = (A 1 ( 2 ;)A 2 (y)A 3 (z)) 1 /3 = fflflfl) '\ 

\Pl P2 P3 J 

(32) 
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III. PHYSICAL AND COMPUTATIONAL PARAMETERS 

The different large-eddy simulations presented here refer to 
the same geometry — see Fig. Q] — and physical parameters as 
the direct numerical simulation (DNS) performed by Leriche 
and Gavrilakis |5j. The details relative to these parameters are 
gathered in Table U The Reynolds number based on the max- 
imum velocity on the lid was chosen to be Re = Uolhjv = 
12 000. 

The kinetic energy is provided to the flow by the shear stress 
at the top lid through viscous diffusion. The amplitude of the 
Reynolds stress below the lid is negligible indicating that the 
flow under the lid is mainly laminar but transient. The mo- 
mentum transfer from the lid induces a region of strong pres- 
sure in the upper corner of the downstream wall as the flow, 
mainly horizontal prior the corner, has to change direction and 
moves vertically downwards. This sharp turn dissipates en- 
ergy in that region. Along the downstream wall the plunging 
flow behaves like a wall jet with a variable thickness. Near the 
symmetry plane the jet thickness is reduced while it increases 
away from this plane. This jet, laminar and unsteady at the 
very beginning, separates from the cavity wall at mid-height 
and grows as two elliptical jets on both sides of the symmetry 
plane. They hit the bottom wall where they produce turbu- 
lence. This turbulence is convected away by the main central 
vortex towards the upstream wall where the flow slows down 
and relaminarizes during the fluid rise. 

Domain size (x,y,z) (2h,2h,2h) 

Wall positions x, y, z = ±h 

Reynolds number Re = U 2h/v 12 000 

No. of spectral elements (E x ,E y ,E z ) (8,8, 8) 

Polynomial orders (N x ,N y ,N z ) (8,8, 8) 

Time step 0.002 h/U 

No. of time iterations 387 000 

Dynamic range 774/t/t/o 

Nodal filtering - DSM & DMM M = N - 2 = 6 



Modal filtering - DSM & DMM (1 st level) k c = N 
Modal filtering - DSM & DMM (2 nd level) k' c = N 



TABLE I: Numerical and physical parameters of the simulations 

In order to resolve the boundary layers along the lid and 
the downstream wall, the spectral elements are unevenly dis- 
tributed as can be seen in Fig. [2] The spatial discretization has 



E 



En 



E z 



3 elements in the three space directions 
with N x = N y — N z = 8 polynomial degree, equivalent 
to 65 3 grid points in total. The spectral element calculation 
has two times less pointsper space direction than the DNS of 
Leriche and Gavrilakis ||5[] who employed a 129 3 Chebyshev 
discretization. Both nodal and modal filters were used in our 
LES computations based on DSM and DMM; the former with 
A I = N — 2 to stabilize the velocity field at each time step 
and the latter with k c = N - 2 (resp. k' c = N — 3) to fil- 
ter the highest modes in the modal Legendre space at the first 
level (resp. second level) of explicit filtering. These filtering 
levels refer to the overbar and the test filtering respectively. It 
is noteworthy recalling here that the modal filter introduced in 
Sec. Ill Dl is not projective. The computations are particularly 
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FIG. 2: Spectral element grid in the mid-plane z/h = 0. 



sensitive to the values of M and k c ; smaller values will affect 
spectral convergence whereas higher values will have very lit- 
tle effect on the smallest scales of the problem. The reference 
results are the DNS data of Leriche 0,151] and the experimental 
ones from Koseff and Street @], corresponding to 1 000 and 
145.5 time units respectively. In the cavity flow, the average 
is obtained by time averaging. 

The LES-DSM and LES -DMM were both started from the 
same initial condition, namely the velocity field obtained from 
the DNS by Leriche and Gavrilakis and re-interpolated from 
the Chebyshev grid onto the spectral-element GLL grid. 

Non-dimensionalization is performed using h as length 
scale, h/Uo as time scale and Uq as velocity scale. All the 
results and data presented in the sequel will be based on this 
non-dimensionalization. 



A. Statistical ensemble averaging 

For any variable, the Reynolds statistical splitting intro- 
duces the average value denoted by a capital letter into brack- 
ets "(•)" whereas a lower case letter will be used to denote 
its fluctuating part. It is noteworthy reminding here the filter 
splitting introduced in Eq. (O using the overbar and prime no- 
tations to denote respectively the resolved and subgrid scales. 
To simplify notations, and unless otherwise stated, the over- 
bar will be omitted in the sequel as most of the fields consid- 
ered are resolved fields derived from solutions of the filtered 
Navier-Stokes equations ([T])-©. More precisely considering 
any variable X can be decomposed as follows 



X= (X) +x = {{X)+x) + {{X')+a/), 



(33) 



where (X) (resp. (X'}) is the average resolved (resp. sub- 
grid) part of X and x (resp. x') is the fluctuating resolved 
(resp. subgrid) part of X. The subgrid scales being unknown, 
the term (X 1 ) + x' cannot be directly computed from the sim- 
ulation. All the results presented in this article refer to re- 
solved quantities be them average (X) or fluctuating x. For 
the sake of simplicity, these quantities are directly and respec- 
tively compared to (X) and x, obtained from reference results, 
see Sec. lIiTCl 
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We assume that a statistically-steady state is attained and 
time averaging will be taken as ensemble averaging. The 
whole dynamic range — cf. Table [Q — corresponding to 1 290 
equally spaced samples has been considered when averaging. 
As the starting point of all LES is the same DNS sample taken 
from a statistically-steady state, it is reasonable to also as- 
sume that these simulations will reach a statistically-steady 
state very quickly. These assumptions can be verified in sev- 
eral number of ways. First, we present in Fig. |3]the time his- 
tories of the volume integral of the total kinetic energy K.(t) 
and the volume integral of the fluctuating energy K,(t) for the 
DNS and both the LES-DSM and LES-DMM. On this figure, 
one can observe that after approximately 80 h/Uo time units, 
the two LES models DSM and DMM start being effective and 
providing different macroscopic results. Both JC(t) and n{t) 
have different time evolutions but within the same range of 
fluctuations and with very close average values, see Table HTl 



UTTl and are showing to be of the order of the error introduced 
by the space and time discretizations. 

Variable Rel. diff. DSM Rel. diff. DMM Anti-/Symmetry 



(U) 


4 


807e- 


04 


6 


696e- 


05 


S 


(V) 


4 


591e- 


04 


3 


014e- 


■04 


s 


(W) 


6 


966e- 


05 


4 


129e- 


04 


A 


(P) 


1 


120e- 


04 


7 


333e- 


05 


S 


s/W) 


8 


758e- 


05 


8 


899e- 


■05 


S 




1 


696e- 


03 


7 


764e- 


04 


S 


vV) 


4 


501e- 


04 


8 


447e- 


05 


S 


(uv) 


1 


107e- 


04 


2 


236e- 


04 


S 



TABLE III: Quantitative assessment of the symmetry and anti- 
symmetry properties of some resolved average fields in the cavity; 
"Rel. diff." stands for maximum relative difference 




100 200 300 400 500 600 700 
non-dimensional time t 

FIG. 3: Time histories of IC(t) (graphs (a)) and of K,(t) (graphs (b)) 
for the DNS (green lines), the LES-DSM (red lines) and the LES- 
DMM (blue lines) 



Average integral terms Magnitude in Uq units 



(tC(t))DSM 
(K.(t))DMM 

(k(£))dns 
(k(£))dsm 
{^(^))dmm 



0.055527 
0.056296 
0.056194 

0.004529 
0.004960 
0.004864 



TABLE II: Average values of K. and k for the DNS, LES-DSM and 
LES-DMM 

A second way to assess the accuracy of the ensemble aver- 
aging is done by testing the property of symmetry (resp. an- 
tisymmetry) with respect to the mid-plane z/h = 0, of some 
first- and second-order statistics of the resolved velocity and 
pressure fields. For each grid point, the relative difference 
between the nodal value at this point and the corresponding 
nodal value at the symmetric grid point is calculated. In the 
antisymmetric case, the opposite nodal value is considered at 
the symmetric grid point. The z-component of the average 
resolved velocity field (W) is the only field presented being 
antisymmetric with respect to the mid-plane z/h = 0. The re- 
sults of the maximum errors on the grid are gathered in Table 



B. Under-resolved DNS and Smagorinsky model 

Before providing the reader with a comprehensive review 
of results obtained for the two models LES-DSM and LES- 
DMM, partial results for the UDNS and the LES-SM are pre- 
sented in this section. These results correspond to the same 
parameters as the one in Table U except that the number 
of iterations is 33 000 corresponding to a simulation length 
of 66 h/Uo — approximately one tenth of the total simulation 
time of the LES-DSM and LES-DMM. Moreover, for the 
LES-SM the value of the Smagorinsky constant Cs defined 
in Eq. © was taken equal to its theoretical value 0.18, see 
lfl2tl for greater details, and no wall-damping procedure was 
implemented for these preliminary simulations. The reference 
result is the DNS by Leriche and Gavrilakis jj] and is rep- 
resented by the solid line in the profiles in Figs. [4] and [5] 
whereas dashed (resp. dotted) lines refer to the UDNS (resp. 
LES-SM). The results in Figs. [4] and E] are one-dimensional 
profiles of the average velocity field and its rms-fluctuations 
in the mid-plane z/h = 0. General conclusions can be drawn 
from all these figures. First, UDNS is totally inoperative in the 
particular context of this simulation. Even first-order statistics 
such as (U) and (V) are far from being well predicted, not 
to mention rms-fluctuations. Second, the Smagorinsky model 
LES-SM results show a real improvement in predicting the 
fields compared to the UDNS but as already known, the sim- 
plicity of this model does not allow to correctly predict the 
stiff physics of this simulation. These results justify the need 
for a more complex LES modeling such as LES-DSM and 
LES-DMM presented in the sequel. 



C. Comparisons with available results 

In this section, results of the LES-DSM and the LES-DMM 
are compared with the available reference experimental and 
numerical results. 
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FIG. 4: In the mid-plane z/h = 0: (U) on the horizontal centerline 
y/h — (top), (V) on the vertical centerline x/h — (bottom): 
DNS (solid line), MILES (dashed line) and LES-SM (dotted line) 



1. One-dimensional profiles 

Of the previous work available in the literature on the lid- 
driven cubical cavity flow, the numerical DNS data from 
Leriche and Gavrilakis Jit], Leriche ll36ll and the experimental 
data of Prasad and Koseff |0] constitute the two main refer- 
ences. This work being an extension of the one by Leriche, it 
borrows from the values of the main physical parameters — 
see Table U The work from Prasad and Koseff ]7[] includes 
data from a flow at Reynolds number similar to that of the 
present LES. The measurements that these authors reported 
were taken in the mid-plane z/h = 0, which is a statistical 
symmetry plane of the flow domain. As it will be shown in 
the sequel, the flow near the downstream secondary eddy — 
see Fig. Q] — is not homogeneous in the z-direction. In the 
"turbulent" part of the cavity, the mid-plane is found to cut 
through surfaces of local minima in the intensity field with 
rapid changes occurring on both sides of it. 

The set of experimental data corresponding to a Reynolds 
number Re = 10 000 is used for the comparisons of the one- 
dimensional average velocity profiles along the vertical and 
horizontal symmetry axes. It is important to note that no ex- 
perimental error-bars were given for any data. The only infor- 
mation related to the local experimental measurement error 
is provided by the two crosses corresponding to two differ- 
ent measurements in the middle (x/h = or y/h = 0) of 
each centerline — the velocity probing system going back and 
forth from this point [13711 . In addition the experimental data of 
Prasad and Koseff iTtL l37h are obtained over a non-dimensional 
averaging time of 145.5 whereas the DNS (resp. LES) results 
were obtained over an averaging time of 1 000 (resp. 774). In 



0.060 




FIG. 5: In the mid-plane z/h — 0: \J (u 2 ) on the horizontal cen- 
terline y/h — (top), yj (v 2 ) on the vertical centerline x/h = 
(bottom): DNS (solid line), MILES (dashed line) and LES-SM (dot- 
ted line) 



absence of local error-bars in the measurements, this may ex- 
plain the scattering (and possible non-convergence) of some 
experimental data, together with practical difficulties of ac- 
curately measuring fluctuating fields in region of low or al- 
most constant velocity. A detailed analysis of the disparity 
between the numerical results and some experimental data can 
be found in 0H. 

For the sets of experimental and DNS data, the total ve- 
locity field is considered whereas in the case of LES, only 
its resolved part is presented. The legend for Figs. [6l4TOl is 
as follows: crosses refer to the experimental data of Prasad 
and Koseff, the solid lines to the DNS by Leriche and Gavri- 
lakis, the dashed lines to the LES-DSM and the dotted lines 
to the LES-DMM. All the data related to average and rms- 
fluctuations of the velocity field are expressed in terms of the 
velocity scale U~o and (uv) in terms of Uq. 

A discussion on the comparisons between the DNS refer- 
ence results and the experimental ones is available in Jit] . In 
the sequel, we will focus on comparing the LES-DSM and 
LES-DMM results with the DNS and experimental ones. A 
rapid overview of Figs. |6l-[T()l indicates that both LES mod- 
els provide results very close to the DNS references, even for 
the rms-fluctuations in Figs. [8] and [9] and above all for the 
component (uv) of the Reynolds stress in Fig. QJJ The dif- 
ferences between the profiles of the two LES models and the 
DNS generally coincide with the existence of local extrema; 
maxima tend to be slightly over-estimated in the LES whereas 
minima are somewhat under-estimated. These two effects can 
be partly justified by the reduced sampling in the LES-DSM 
and LES-DMM compared to the sampling of the DNS. This 
phenomenon have already been encountered and studied by 
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FIG. 6: In the mid-plane z/h = 0: (U) on the horizontal centerline 
y/h — (top), (U) on the vertical centerline x/h — (bottom); Ex- 
periment (crosses), DNS (solid line), LES-DSM (dashed line), LES- 
DMM (dotted line) 
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FIG. 8: In the mid-plane z/h = 0: \J (u 2 ) on the horizontal center- 
line y/h = (top), y ' (u 2 ) on the vertical centerline x/h — (bot- 
tom); Experiment (crosses), DNS (solid line), LES-DSM (dashed 
line), LES-DMM (dotted line) 
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FIG. 7: In the mid-plane z/h = 0: (V) on the horizontal centerline 
y/h — (top), (V) on the vertical centerline x/h — (bottom); Ex- 
periment (crosses), DNS (solid line), LES-DSM (dashed line), LES- 
DMM (dotted line) 
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FIG. 9: In the mid-plane z/h = 0: ^/ (v 2 ) on the horizontal center- 
line y/h — (top), \J (v 2 ) on the vertical centerline x/h = (bot- 
tom); Experiment (crosses), DNS (solid line), LES-DSM (dashed 
line), LES-DMM (dotted line) 



Leriche in 0]. From these results it is not possible to rank be- 
tween themselves the performances of the LES-DSM and the 
LES-DMM. 



2. Two-dimensional profiles 

The comparisons with the DNS results started in Sec. 
IIII C 1 1 are now extended to the whole mid-plane z/h = by 
plotting identical series of contour levels of average velocity 
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FIG. 10: In the mid-plane z/h — 0: (uv) on the horizontal center- 
line y/h — (top), (uv) on the vertical centerline x/h — (bot- 
tom); Experiment (crosses), DNS (solid line), LES-DSM (dashed 
line), LES-DMM (dotted line) 



components in Fig. [TTJand of rms-fluctuations of the velocity 
in Fig. El for the DNS (left column), the LES-DSM (central 
column) and the LES-DMM (right column). 

As previously noted with the one-dimensional profiles, the 
results provided by the LES-DSM and LES-DMM are both 
very close to the reference DNS results. Secondary corner 
eddies located above the bottom wall and below the lid next 
to the upstream wall are correctly captured in the mean flow. 
Other finer structures visible in Fig. [T2] (bottom), for \J (v 2 ) 
near the upstream wall are also correctly captured by both 
LES modelings. The rms-fluctuations of the x-component u 
of the velocity field is accurately resolved just below the lid 
which is a high-gradient region for the mean flow. Moreover, 
in the region near the downstream wall where the wall jet — 
separated into two elliptical jets — are impinging on the bot- 
tom wall, the high gradients of velocity fluctuations are well 
reproduced. As it will be presented in the following sections, 
the maximum of turbulence production belongs to this region 
of the flow domain which will be indeed of particular interest 
in the remaining. 

The flow below the lid and near the corner with the down- 
stream wall presents wiggles in the LES contours for (V), see 
Fig. QT| (bottom). Although less intense, these wiggles are 
also noticeable on the contours for J (v 2 ) at the same loca- 
tion, see Fig. [12] (bottom). More limited effects are noticeable 
for the equivalent .T-component fields. These very limited de- 
fects in both simulations find their origin in a slight under- 
resolution of the spectral-element grid in this small region of 
the cavity where high gradients are present. 



D. Physical parameters of the LES modeling 



The LES modeling for both the LES-DSM and LES-DMM 
involves the calculation of two scalar fields, namely the dy- 
namic parameter Ci and the eddy viscosity v sgs which are 
inter-related. As some values of Cd produced by the two dy- 
namic procedures (0 and (Tl2l may locally reach relatively 
"high values" destabilizing the time-integration procedure for 
the filtered Navier-Stokes computations. Hence it is common 
to use ad hoc averaging or limiting of the dynamic parame- 
ter Cd to ensure stability. Various procedures are reported in 
the literature: averaging in homogeneous directions lfl4ll38ll . 
temporal smoothing 113911 . integral constraint l40ll . Lagrangian 
averaging BUl and clipping iflil l29ll . In the present work the 
latter procedure of clipping is used. First the maximum ad- 
mitted value of the dynamic parameter was Cd max = (0.18) 2 , 
0.18 corresponding to the theoretical value of the Smagorin- 
sky constant — see lfl2ll . The negative values of Cd are also 
clipped and set to zero for the LES-DSM and LES-DMM. The 
amount of grid points clipped is indeed very limited and cor- 
respond to 0.2% and 0.08% of the total number of grid points 
for LES-DSM and LES-DMM respectively. It was found that 
the clipping of Cd to the interval comprised between — (0.18) 2 
and +(0.18) 2 — therefore allowing for local negative values of 
the eddy viscosity — was not affecting at all the stability of the 
spectral-element filtered Navier-Stokes computation. The dif- 
ference between the results with or without the negative val- 
ues of Cd was found to be negligible in the particular context 
of the lid-driven cubical cavity flow which is related to the 
limited amount of backscattering for this flow at a Reynolds 
number of 12 000. 

Fig. [T3]displays contour lines of the average eddy viscosity 
for the LES-DSM in the mid-plane z/h = and in the plane 
z/h = 0.241 where the maximum of average turbulent en- 
ergy dissipation rate was localized — cf. Sec. II V CI for greater 
details. First, the C°-continuity breakage in the inter-element 
continuity is obvious — see Fig. [2]to compare with the spectral 
element grid in the mid-plane — and is directly related to the 
discontinuous nature of the filter length field A defined in Sec. 
Ill D 21 Eq. d32l . The effect of such discontinuity of the sub- 
grid viscosity has been analyzed and discussed by Blackburn 
and Schmidt in j29ll using the same numerical framework as 
ours, namely the SEM. They found that the inter-element dis- 
continuity of the subgrid term does not have a noticeable ef- 
fect on their physical results which is confirmed by the present 
work. Finally, it appears clearly that the reasons for resorting 
to a dynamic procedure are fully justified by Fig. [13] Indeed, 
the dynamic procedure automatically turns on the dynamic pa- 
rameter Cd which in turn activates subgrid-scale viscous ef- 
fects in the regions of the flow where turbulent dissipation at 
the small-scales level occurs — see Sec. |IV C| 

Similar results are obtained for the eddy viscosity and the 
dynamic parameter in the case of the dynamic mixed model 
LES-DMM. The same clipping procedure, with the same clip- 
ping values as described earlier for LES-DSM was imple- 
mented. 
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FIG. 12: Contours of rms-fluctuations of the velocity in the mid-plane z/h = 0; DNS (left), LES-DSM (center) and LES-DMM (right); 20 
contours equally spaced between and 0.1 for \J (u 2 ) (top) and between and 0.15 for -J (v 2 ) (bottom) 
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FIG. 13: Contours of the average eddy viscosity {v sgs ) for the LES- 
DSM in the mid-plane z/h = (left) and in the plane z/h = 0.241 
(right); same series of contour levels is used in the two planes 



IV. CHARACTERIZATION OF TURBULENCE IN THE 
FLOW 



which in turn can be recast in terms of u> the fluctuating vor- 
ticity, fi being the total resolved vorticity field 



(e) = v{LOiU)i) + 2u 



d 2 (ujUj) 
dxidxj 



(36) 



We define the average difference (S) by the difference be- 
tween the average turbulent energy dissipation rate, divided 
by v, and the average fluctuating enstrophy 



(uJiUJi) 



dxidx-j 



(37) 




This section is devoted to a thorough analysis of some spe- 
cific features of the flow in the region of the cavity where tur- 
bulence occurs. The aims are to ensure that the LES-DSM and 
LES-DMM are both capable of reproducing the fine physics 
observed in these regions and also to gain insights in the tur- 
bulent mechanisms involved. 



A. Inhomogeneity of turbulence 

It is easily predictable that such a confined flow will pro- 
duce an inhomogeneous turbulence but it is worth determin- 
ing in greater details the turbulent inhomogeneous zones in the 
cavity. In order to access this information we use the average 
turbulent energy dissipation rate (e) defined by 



/ dui 


h du A 


( dui 




\dxj 


dxi ) 


\dx 3 





2v(SijSij) 



(34) 

Here and in the sequel, we use index notation and the sum- 
mation convention, where repeated indices imply summa- 
tion. The velocity fluctuations being divergence-free, one can 
rewrite 



dui dui 
dxj dxj 



d 2 (ujUj) 
dx;dxj 



(35) 



FIG. 14: Region of the cavity where the turbulent flow is inhomoge- 
neous according to the criterion \{8) |/|{5)max| > 1/100; LES-DMM 



For homogeneous flow, the spatial derivatives of the 
Reynolds stress components (iMUj) are zero, and subse- 
quently (S) = 0. The average difference (6) was calculated 
for both databases LES-DSM and LES-DMM. Fig. Qj] dis- 
plays a 3D view of the volume of the cavity where the flow is 
inhomogeneous according to the following heuristic criterion: 
IvWIWmaxI > 1/100, where \{5)™*\ = 159.8U$/h 2 and 
IW™ M | = 155 - 6 U$/h 2 . In other words, it shows the re- 
gion of the flow where the inhomogeneity of the turbulence — 
measured by (S) — is above 1% of its maximum absolute 
value. 

As expected, one can observe in Fig. [15] that in the region 
near the downstream wall where the two primary elliptical jets 
are impinging on the bottom wall, the flow is highly inhomo- 
geneous. More specifically, the inhomogeneity is more im- 
portant in the zone in between the two elliptical jets where the 
flow is ejected and recirculating. Likewise similar patterns 
with lower magnitudes are detected in the regions where the 
secondary jets and the tertiary jets are impinging. The sec- 
ondary jets are impinging on the bottom of the upstream wall 
producing an inhomogeneous turbulence visible in Fig. [T5lfor 
values of x/h close to — 1 . For the tertiary jets, impinging on 
the upstream part of the lid, the inhomogeneity is only visible 
in the 3D view in Fig. [14] 



13 




-S35S : 

'-1 -0.5 „ 0.5 1 -0.5 „ 0.5 1 



FIG. 15: Contours of (8) in the plane y/h = —0.968 just above 
the bottom wall (left) and in the mid-plane z/h = (right); 100 
equally spaced contours corresponding to levels between the thresh- 
old 0.01|(5) ra ax|; dashed contours correspond to negative levels with 
a colormap ranging from blue to red; LES-DMM 

B. The turbulence production near the downstream wall 

As mentioned by Leriche and Gavrilakis in the largest 
turbulence production rates in the cavity are to be found in the 
primary elliptical jets parallel to the downstream wall, near the 
impact points just above the bottom wall. The budget equa- 
tions of the resolved second-order moments (uiUj) govern- 
ing the turbulence energetics — see l42tl and l43ll for greater 
details — comprise a term named here Py , defined by 

Pij = . M ?m. M m., (38) 

and corresponding to the interaction of the resolved mean flow 
and the resolved Reynolds stress tensor. P^ can be interpreted 
as responsible for the production of Reynolds stresses or in 
other words for the production of turbulence. 

1. Maximum of turbulence production near the downstream wall 
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FIG. 16: Contours of resolved the production of turbulence term 
P22 in the plane y/h = -0.9388; LES-DSM (top) and LES-DMM 
(bottom); 20 contour levels equally spaced between —0.025 Uq /h 
and 0.070 U^/h; dashed lines refer to negative contour levels; bullet 
points • refer to the same grid node of coordinates (x/h = 0.7874, 
y/h = -0.9388, z/h = 0.3371) 



In the specific case of the separated downstream-wall jet, 
the term P22 is the largest out of the set of turbulence pro- 
duction terms {Pij}. After probing in the cavity, the max- 
ima of the resolved field P22 is to be found in the plane 
y/h = — 0.9388just at a very short distance above the bottom 
wall. The maximum values obtained are P^ SM = 0.070 Uq /h 
and P 2 D 2 ^ = 0.064 U$/h. The contours of the turbulence 
production term P22 in this plane y/h = —0.9388 are shown 
in Fig. [16] for both LES models. First, it can be noted that 
these contours are qualitatively very close to the ones obtained 
by Leriche and Gavrilakis in Jit]. For x/h > 0.5, the distri- 
bution of contours of the production of turbulence allow to 
clearly visualize the trace of the separated elliptical jets. 

2. Time histories and power spectra at the maximum of turbulence 
production 

In Fig. [16] one can notice that for each LES model, one 
grid point — identical for both DSM and DMM — has been 



highlighted with a bullet point •. This point denoted by Go 
whose coordinates are x/h = 0.7874, y/h = —0.9388, 
z/h = 0.3371, is the closest grid point to the two maxima 
of P 22 for LES-DSM and LES-DMM. The point 6 provides 
the optimal search position for probing time histories of vari- 
ous turbulent fields in the sequel. 

First, the values of the 2>component of the fluctuating re- 
solved velocity field u, of the fluctuating resolved pressure p 
and the resolved turbulent kinetic energy k = mUi/2, have 
been extracted of the LES-DSM database for each and every 
sample. These time histories are shown in Fig. [17] — note that 
only the last 1 024 samples out of the total of 1 290 that con- 
stitutes the database are presented. 

Based on these results the corresponding power spectra 
have been computed by fast Fourier transform — a posteriori 
justifying the choice of 1 024 samples in the previous time 
histories — and are presented in Fig. [18] The scattering of 
points in the high-wavenumber Hi zone is expected for spectra 
of non-spatially average fields. For such inhomogeneous flow 
with highly localized turbulent effects, averaging in space the 
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FIG. 17: Time histories of p (graph (a) shifted of +0.25), the turbu- 
lent kinetic energy k (graph (6)) and it (graph (c) shifted of —0.25), 
at the point of coordinates (0.7874, -0.9388, 0.3371); LES-DSM 
database 
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FIG. 18: Power spectra for the fluctuating pressure p (a), the tur- 
bulent kinetic energy k (b), and the fluctuating velocity field u (c), 
obtained from the time histories in Fig. [T7J LES-DSM database 



fields not only pleasantly reduce the scattering of points but 
concurrently strongly modifies the high-wavenumber scaling 
which is the main source of information brought by the spec- 
tra. Nevertheless, the spectra offer a qualitative information 
regarding the Eulerian time scales of the spatial structures of 
turbulence convected past the point 8o- The resolved mean 
flow depicted in Fig. [TJ] — highlighting the presence of the 
core central primary vortex and secondary corner vortices — 
serves merely to convect turbulence at the bullet-spotted point 
0o. A careful scrutinizing of the resolved mean flow in the 
vicinity of 0o shows that this point is exactly positioned at 
the "interface" between the core central vortex and the bottom 
corner vortex. The power spectra in Fig. Q~8]feature the distri- 
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FIG. 19: Two-dimensional projected average resolved velocity vec- 
tors in the mid-plane z/h = (top) and in the plane y/h = —0.9388 
(bottom); bullet point refers to Go; LES-DSM; bullet points • refer to 
the same grid node of coordinates (x/h = 0.7874, y/h — —0.9388, 
z/h = 0.3371) 

bution of frequencies — or equivalently of time scales. These 
frequencies ui refer to the convection past 0o of turbulent 
structures of size of order I at a velocity of order (U)(@o), 
leading to the relation 



(u)(e ) 



(39) 



The average resolved velocity field at 0q being given, the 
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spectra hence instruct us on the distribution of spatial scales 
of resolved turbulent structures convected by the mean flow 
at this point where the production of turbulence is maximum. 
Unfortunately, the relatively low sampling of the LES-DSM 
database and the not-long-enough simulation range interval 
does not permit to reach the highest frequencies of the order_, 
of (f7)(0 o )/A where A is the filter length— see Eq. d32|i— 
defining the LES scale separation. 



3. Determination of coherent structures responsible for the peaks 
of turbulence production 

In an attempt to provide a comprehensive and thorough as- 
sessment of the performances of both LES models, the deter- 
mination of the coherent structures responsible for the intense 
turbulence production at the point O has been envisaged as 
an ultimate challenge for both SGS modeling. The first step 
towards this goal necessitates to study the instantaneous dis- 
tribution of the resolved term —v 2 d(V)/dy which was found 
to be the predominant term in P22, see Fig. l20l displays 
the time histories of this term for the DNS, LES-DSM and 
LES-DMM. Both LES present a limited number of high-value 
peaks which are assumed to be engendered by specific coher- 
ent vortices or large eddies. The intensity of the peaks pro- 
duced by the LES-DMM is lower than those generated by the 
LES-DSM. This is supposed to be due to an over-evaluation of 
the eddy viscosity by the dynamic procedure of the DMM. In 
addition, this is consistent with the observation made in Sec. 
II VB II where P™™ < P 2 D ™ was found and with the values 
of resolved (K.) and (n) in Table [TT1 

In order to finally characterize possible large eddies which 
would be responsible for these peaks, database samples pro- 
ducing a resolved term —v 2 d{V)/dy above the threshold 
value 0.15 U^/h were put aside to form a subset of the com- 
plete databases. The size of the subset of samples for LES- 
DSM (resp. LES-DMM) is approximately 6% (resp. 5%) of 
the size of the complete database. Based on these two subsets 
a conditional averaging of the streamwise resolved vorticity 
field fl x is performed. Fig. |2U displays the contours of this 
quantity in the vicinity of Qq — the domain represented corre- 
sponds to only 4% of the surface of a normal section of the 
cavity — for both models. Two counter-rotating vortices are 
clearly exhibited by both models, together with the intense in- 
fluenced shear layers laying on the bottom wall at y/h = — 1. 
This vortex pair is identified as the coherent eddy responsible 
for the turbulence peaks and production in this region. The 
characteristic length scale of this large eddy is of the order 
of O.lh. Having identified this vortex pair we can further an- 
alyze the time histories at 6n of the resolved pressure and 
the resolved turbulent kinetic energy depicted on Fig. [TTl 
graph (a) and (b) respectively. One can notice that the in- 
tense peaks of the resolved term —v 2 d(V)/dy on Fig. |20] 
correspond to intense peaks of resolved turbulent kinetic en- 
ergy on Fig. [T7] (6) and to low-pressure peaks on Fig. [T7] 
(a). The vortex pairs generated by this turbulent flow are 
responsible for the low pressures and the high turbulent ki- 
netic energy thereby justifying the observed correlations be- 
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FIG. 20: Time histories of the resolved term —v 2 d(V)/dy in Uo /h 
units for the DNS (green), LES-DSM (red) and the LES-DMM 
(blue); the dotted lines represent the threshold value 0.15 Uq /h 



tween these three time histories. Moreover the intensity of the 
vortex pair calculated by the LES-DSM is again higher than 
the one from the LES-DMM: (ft° SM ) C a milx = 18.6 U /h and 
(f2° MM ) camax = 14.0 Uo/h, where the subscript "ca" stands for 
conditionally averaged. The intensity, the more regular struc- 
ture and the localization of the vortex pair are three features 
suggesting that the dynamic Smagorinsky model provides a 
better SGS modeling than the dynamic mixed model. Never- 
theless, the DMM performances in terms of SGS modeling are 
more than satisfactory. When considering the complete aver- 
aging of the ^-component of the resolved vorticity field in the 
region where the vortex pair has been localized by conditional 
averaging, see Fig. [2T| (n x ) was found almost constant and 
of magnitude approximately 0.9 Uo/h. 
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C. Small-scales turbulent structures 



A characteristic of high-Reynolds-number turbulence is 
that the vorticity possesses intense small-scale, random vari- 
ations in both space and time. The spatial scale for vortic- 
ity fluctuations is the smallest in the continuum of turbulent 
scales, i.e the Kolmogorov scale. Analogously to the vortic- 
ity fluctuations, for large-Reynolds-number turbulence veloc- 
ity gradients dui/dxj are also dominated by the small scales 
of turbulence and the overall energy dissipation rate of kinetic 
energy is dominated by the average turbulent energy dissipa- 
tion rate (e) defined in Eq. d34i i. In the LES framework, the 
interest for small scales is twofold. First, small scales fall into 
the range of subgrid scales and therefore are not simulated but 
wholly modeled to properly reproduce their interactions with 
larger resolved scales of the flow. Second, the small scales 
have the crucial role to terminate the turbulent energy cascade 
by dissipating the energy originating from large eddies. An in- 
correct SGS modeling will produce either an over-dissipation 
or conversely an under-dissipation of kinetic energy. The time 
histories of the total kinetic energy of the cavity flow IC(t) 
and of the total turbulent fluctuating energy K(t) presented in 
Fig. [3] are in this framework a precious proof of the correct 
global prediction of the energy dissipation by the modeled 
small scales in volume. 



1. Localization of small-scales structures 



FIG. 22: Visualization of the region of the cavity where the average 
turbulent energy dissipation rate (e) is above 1% of its maximum 
value 3 570 vUljb?; LES-DMM 



is worth keeping in mind that the more intense (e) the more 
small scales are involved in the dissipation process. Fig. [23] 
displays with decreasing intensity, the dissipation due to the 
impingements of the separated wall jets on the bottom wall 
(bottom-left), on the upstream wall (bottom-right) and on the 
lid-plane (top-left). It appears that the LES-DSM is not able 
to properly reproduce the same intensity for the two symmet- 
ric jets impinging on the upstream wall (bottom-right). The 
same asymmetry in the intensity of (e) is observed for the 
LES-DMM which could presumably be due to the observed 
asymmetry — with respect to the mid-plane — of the eddy vis- 
cosities generated by the dynamic procedures of both SGS 
modeling. 

In Fig. [23] (bottom-left), one can notice that one grid point 
has been highlighted with a black rectangle □. This point 
denoted by So whose coordinates are x/h = 0.7685, y/h = 
—l,z/h= 0.2410, is the closest grid point to the maximum of 
(e) for LES-DSM. The point So provides the optimal search 
position for probing small-scales related fields. The plane-cut 
z/h = 0.241 — passing by So — of (e) in Fig. [23] (top-right) 
exhibits a qualitative correlation with the same plane-cut for 
the average eddy viscosity (f sgs } in Fig. [13] (right). 



In this context, it appears relevant to first locate small-scales 
turbulent structures in the cavity and afterwards to check the 
correlation between the small-scales positioning and the ac- 
tivation of the SGS modeling represented in Fig. [13] Small 
scales can be indirectly localized by investigating the zones 
of intense average turbulent energy dissipation rate. Indeed 
(e) involves products of fluctuating velocity gradients, see Eq. 
fl34l l. First qualitatively, the region of the cavity flow corre- 
sponding to values of (e) above 1% of its maximum value is 
shown in Fig. [22]for the LES-DMM. As foreseen, the wall-jet- 
impinging regions are subject to intense turbulent energy dis- 
sipation at the small-scales level. The two-dimensional cuts 
in Fig. |23l offer a more detailed information regarding the in- 
tensity of (e) in four different planes of specific interest. It 



2. Correlation between small-scales localization and eddy 
viscosity 

Such correlation between the small-scales localization and 
the activation of the LES dynamic Smagorinsky modeling is 
important in practice to ensure the effectiveness of the SGS 
modeling. Therefore a more quantitative approach is required, 
which relies on the calculation of a correlation field based on 
the instantaneous values of e and is sgs . The following correla- 
tion coefficient C, defined by 



C = C(e, K;g S ) 



(e)<*V) 



«£" 2 ><^» 1/2 



(40) 
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FIG. 23: Two-dimensional contour lines of (e) in the following 
planes: lid-plane y/h — 1 (top-left), plane z/h = 0.241 (top-right), 
bottom-plane y/h = —1 (bottom-left), upstream-plane x/h = —1 
(bottom-right); LES-DSM; black rectangle □ refers to the grid node 
of coordinates (x/h = 0.7685, y/h = -1, z/h = 0.2410) 



where " stands for the fluctuating part of the considered field, 
was calculated for the complete set of samples LES-DSM. 



calization. More precisely the poor correlation in the vicinity 
of S is imputed to the fact that in this region, the term Sy 
in the turbulent energy dissipation rate — see Eq. (f34-b — varies 
very rapidly in space likewise the eddy viscosity. At this point, 
the information provided by the analysis of the subgrid-scale 
activity in the next section is a good complement to the previ- 
ous correlation study. 



3. Subgrid-scale activity 

The filtered kinetic energy can be decomposed into the ki- 
netic energy of the resolved velocity field and the residual ki- 
netic energy which is equal to The conservation equa- 
tion for the kinetic energy of the resolved velocity field — see 
pp. 585-586 in ll43ll for greater details — comprises transport 
terms as well as source/sink terms which are of prime inter- 
est. First is the sink term e y — 2vSijSij = 2SijSij/Re 
corresponding to the viscous dissipation associated with the 
resolved velocity field. The second sink term e sgs = —TfjSij 
corresponds to the SGS contribution and represents the rate of 
transfer of energy from the resolved scales of the flow to the 
subgrid scales. This term e ses is often inappropriately referred 
to as the SGS dissipation in the literature. Indeed e sgs does 
not correspond to any physical dissipation but finds its origin 
in inertial processes. In addition, it is important to note that 
locally e sgs can take negative values. 

The SGS activity, denoted by „4 sgs in the sequel, allows to 
study the local energy fluxes due to the SGS effects. Follow- 
ing Geurts and Frohlich in |44| and Meyers et al. in l45ll . _4 sgs 
is defined as 




FIG. 24: Contours of the correlation coefficient C in the plane z/h — 
0.241 containing the point H ; LES-DSM 

Contours of C in the plane z/h = 0.241, passing by So are 
presented in Fig. [24] The high-correlation zones reproduce 
in essence the turbulent-dominated regions of the cavity and 
even suggest the mean-flow convective effect of the central 
core vortex and other secondary corner vortices on the tur- 
bulent pockets. Nevertheless, higher correlations would have 
been expected in the vicinity of So. Such low correlations 
are evidences of the limitations of the LES in this region of 
the cavity flow. Conversely, the high correlations near the up- 
stream wall are in good agreement with the small-scales lo- 



~7~fjSij + 2uSij Sij 



(41) 



The SGS activity ,4 sgs measures the importance of the subgrid 
scales in the overall dissipation process of the kinetic energy 
of the resolved velocity field. As mentioned by Meyers et al. 
in ll45ll . the SGS activity varies between zero and one where 
a value of zero corresponds to DNS and As gs = 1 is asso- 
ciated with LES at infinite Reynolds number. Moreover the 
value of yl sgs is directly related to the filter width A and mea- 
sures the "distance" between a DNS resolving all flow fea- 
tures at sufficiently high spatial resolution and an actual LES 
corresponding to a specific filter width and mesh spacing. In 
the particular case of the LES-DSM, the SGS sink term is 

e SBS = 2v sgs SijSij = 2C d A |S| SySy, leading to 



jDSM _ 
-^sgs — 



C d A 2 \S\/v 
+ C d A 2 |S|/i/ 



(42) 



Fig. [25] displays the average value (.Asgs) of the SGS ac- 
tivity in the plane z/h = 0.241 containing So and where the 
turbulent energy dissipation rate is maximum. First, it appears 
that the SGS activity is slightly higher for the LES-DMM than 
for the LES-DSM. Moreover, it appears very clearly that the 
SGS modeling is activated in the region of the cavity where 
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the different wall jets are present, with maxima in the im- 
pingement zones. The LES-DMM is more effective in acti- 
vating the subgrid scales in these particular zones. In the zone 
where the tertiary wall jet is impinging on the lid, SGS dissi- 
pation for the LES-DSM is less than 25% of the total dissipa- 
tion, whereas it is above 45% for the LES-DMM. 
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FIG. 25: Contours of the average SGS activity (-4 sgs ) in the plane 
z/h = 0.241 containing the point H : LES-DSM (top) and LES- 
DMM(bottom); same series of contour levels is used for both models 
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FIG. 26: Contours of the average resolved helicity density (h) in 
the bottom plane y/h = —1 (left) and in the plane x/h = 0.7874 
containing Bo; LES-DSM 



Mappings of the average resolved helicity density in Fig. [26] 
allows to locate resolved helical coherent structures (HCS). 
These HCS are particularly intense in the secondary-corner- 
eddy region and are consistent with the experimentally ob- 
served typical HCS, namely streamwise counter-rotating vor- 
tices 14711 . This pairing of coherent helical structures corre- 
spond to a pairing of coherent vortical structures having oppo- 
site vorticity and consequently opposite helicity. Such obser- 
vation justifies the relatively small — but non-zero — resolved 
average helicity reached by both LES models: (H DSM ) = 
0.00764 U$h 2 and (H DMM ) = -0.00 572 U$h 2 . Smaller HCS 
have been identified earlier in Sec. IIVB 31 where stream- 
wise counter-rotating vortices — cf. Fig. [2J] — near the bot- 
tom wall, have been identified by the conditional averaging 
as the principal coherent structures responsible for the high- 
intensity peaks in the production of turbulence in this region 
of the flow. Finally, it is noteworthy to emphasize the strong 
link between average resolved helicity density contours and 
average resolved turbulent kinetic energy dissipation rate ones 
in Fig. [23] 
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D. Helical properties of the cavity flow 

The helicity Ti of the fluid flow confined in the cavity V at 
instant t is defined by 



H(t) 



u ■ uj dV, 



(43) 



and is a measure of linkages and knots between the vorticity 
lines of the flow. The quantity h(x, t) = u • u> is the helicity 
density and is a pseudo-scalar quantity just like TL. The he- 
licity is an important flow quantity because just like the total 
energy of the flow /C, it is an invariant of three-dimensional 
homogeneous turbulence J46ll . The study of the resolved he- 
licity TL and the average resolved helicity density (h) in the 
particular context of the lid-driven cavity flow in a locally- 
turbulent regime allows to gain insights into very important 
features of the turbulence dynamics ll46Tl . For instance TGL 
vortices and secondary corner eddies are structures encoun- 
tered in the lid-driven cavity flow which are well known as 
typical helical structures. 
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non-dimensional time t 

FIG. 27: Time histories of H(t) for the LES-DSM (red lines) and the 
LES-DMM (blue lines) 

As mentioned previously the average resolved helicity of 
the flow is non-zero. The time histories of the resolved helic- 
ity for the whole simulations are shown in Fig. [27] for both 
LES-DSM and LES-DMM. In Sec. ImAl was mentioned that 
both SGS models start being effective and producing differ- 
ent global results after a transient period of about 80h/Uo 
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time units, and likewise the helicity as can be seen in Fig. 
[271 Moreover, the amplitude of the resolved helicity fluctua- 
tions is not decaying during the simulation and the LES-DMM 
qualitatively produces more high-amplitude negative helical 
values therefore justifying its negative average value. 
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FIG. 28: One-dimensional relative helicity spectrum a(k); LES- 
DMM 

Helicity, like energy, is cascaded from large scales down to 
the Kolmogorov dissipation scale, where it is destroyed. Un- 
fortunately, the relatively low Reynolds number of both LES 
does not permit the determination of quantitative scalings of 
energy and helicity spectra which could be compared to the 
Kolmogorov scalings in k~ 5/>3 in helical three-dimensional 
homogeneous isotropic turbulence, as mentioned by Borue 
and Orszag in l48ll . In the same paper, Borue and Orszag con- 
clude that helicity is inherently a large-scale quantity which 
behaves similarly to a passive scalar. Consequently the one- 
dimensional relative helicity spectrum defined by 



a{k) 



H(k) 
2kZ(k) 



(44) 



where H (resp. /C) is the one-dimensional resolved helicity 
(resp. energy) spectrum, decreases at small scales. Even if in 
our context, the turbulence is not homogeneous nor isotropic, 
the previous assertion is undeniably verified by both LES as 
can be seen only for the LES-DMM in Fig. [28] for high val- 
ues of k corresponding to small scales. Similar relative he- 
licity spectrum is obtained for the LES-DSM. This suggests 
that the decreasing trend at small scales of the relative helicity 
spectrum is more general and not only limited to the homo- 
geneous and isotropic turbulence theoretical framework, just 
like the Kolmogorov scale in the inertial range. 



These simulations were based on an accurate spectral-element 
spatial discretization, having two times less points per space 
direction than the direct numerical simulation reference re- 
sult from Leriche and Gavrilakis fH. All filtering levels intro- 
duced in both SGS modelings rely on explicit modal filters in 
the spectral space, retaining C°-continuity of the numerical so- 
lution of the filtered Navier-Stokes equations. An additional 
nodal filter was used to stabilize both LES. Time-averaging 
was shown to be equivalent to ensemble-averaging, with re- 
spect to the global precision level of the numerical integration. 

Partial simulation results using the UDNS and the 
Smagorinsky model as subgrid-scale models, have served to 
prove the necessity of a dynamic SGS procedure. Full LES 
results for both dynamic models have shown very good agree- 
ment with the DNS reference results. The agreement with the 
experimental reference results from Prasad and Koseff [Q] is 
qualitatively good. 

At a Reynolds number of 12 000, the lid-driven cavity flow 
is placed in a locally turbulent regime and such turbulent flow 
is proved to be highly inhomogeneous in the secondary-corner 
regions of the cavity where turbulence production and dissi- 
pation are important. The maximum production of turbulence 
was found to be located in the downstream-corner-eddy re- 
gion just above the bottom wall. An analysis of the spectra of 
turbulent quantities at this point allowed us to determine the 
distribution of the scales of the turbulent structures convected 
past this maximum. Moreover, both LES were able to capture 
the coherent counter-rotating pair of vortices which are mainly 
responsible for the peaks of turbulence production still at this 
point. LES-DSM have shown globally more intense and better 
results than the LES-DMM in this matter. 

Small-scales turbulent structures were located indirectly by 
studying the regions of intense turbulent energy dissipation 
rate e. The eddy-viscosity field was shown to be strongly cor- 
related to e in the turbulent areas of the flow, but the clipping 
procedure — necessary for stabilizing the numerics — imposed 
to the dynamic parameters strongly diminishes this correla- 
tion in the intense turbulent zones. Subgrid-scales activity has 
been analyzed and the higher SGS activity is associated with 
the LES-DMM. 

Helical properties of the flow were investigated. Typical 
helical coherent structures were identified in the secondary- 
corner regions. These structures appear to be correlated to 
the turbulent energy dissipation rate e. The relative helicity 
spectra is shown to be decreasing at small scales, which is in 
agreement with the theoretical results from Borue and Orszag 
14811 for the three-dimensional isotropic inhomogeneous tur- 
bulence. 
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